; Extract several variables at diff. press levels and transfer them to daily scale

; Path
inpath   = "../../hourly_to_daily/"
outpath  = "./"

; Define time
TIME = yyyymmdd_time(2018, 2021, "integer")
time = TIME({20181001:20181231})    ; coordinate subscripting

do d=0,dimsizes(time)-1

t 			= 	sprinti("%0.8i",time(d))

t_name_in = str_get_cols(t,0,3) + "-" + str_get_cols(t,4,5)\
		 + "-" + str_get_cols(t,6,7) + "_" + sprinti("%0.2i",0)\
		 + "_00_00.nc"

print(t_name_in)

fin = addfile(inpath + "wrfout_daily_d02_" + t_name_in, "r")

; ext data
Times        = fin->Times
lat2d		 = fin->XLAT(0,:,:)
lon2d		 = fin->XLONG(0,:,:)
QVAPOR_eta   = wrf_user_getvar(fin,"QVAPOR",-1)     ; The variable to interpolate
wa_eta       = wrf_user_getvar(fin,"wa",-1)
uvmet_eta    = wrf_user_getvar(fin,"uvmet",-1)
height_eta   = wrf_user_getvar(fin,"height",-1)
pw			 = wrf_user_getvar(fin,"pw",-1)
ps    		 = fin->PSFC
ps			 = ps/100 ;hpa

; from uvmet to U and V
uamet_eta   = uvmet_eta(0,:,:,:,:)
vamet_eta   = uvmet_eta(1,:,:,:,:)

; inter from eta level to press level (xy-pl)
vert_coord          = "pres"
interp_levels       = ispan(1000, 100, 20)		; the inter data is consistent with this order
interp_levels@units =  "hpa"
opts             = True
opts@extrapolate = True 		; under ground
opts@field_type  = "pres"
opts@logP        = True   ; Set true when inter press

Q_PL      = wrf_user_vert_interp(fin,QVAPOR_eta,vert_coord,interp_levels,opts)
W_PL      = wrf_user_vert_interp(fin,wa_eta,vert_coord,interp_levels,opts)
Umet_PL   = wrf_user_vert_interp(fin,uamet_eta,vert_coord,interp_levels,opts)
Vmet_PL   = wrf_user_vert_interp(fin,vamet_eta,vert_coord,interp_levels,opts)
height_PL = wrf_user_vert_interp(fin,height_eta,vert_coord,interp_levels,opts)
; delete
delete(uamet_eta)
delete(vamet_eta)

; to specific humidity
Q_PL = Q_PL/(Q_PL+1)
Q_PL@units = "kg/kg"
Q_PL@description = "specific humidity"

; UQ VQ
UQ_PL = Q_PL * Umet_PL
VQ_PL = Q_PL * Vmet_PL
UQ_PL@units = "["+Q_PL@units+"]["+Umet_PL@units+"]"
VQ_PL@units = "["+Q_PL@units+"]["+Vmet_PL@units+"]"
copy_VarCoords(Umet_PL, UQ_PL)
copy_VarCoords(Vmet_PL, VQ_PL)

; Cal vibeta UQ VQ
linlog = 1
g = 9.8
pbot = max(interp_levels)
ptop = min(interp_levels)
UQ_vint = 1./g * 100 * vibeta (interp_levels,UQ_PL(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,ps,pbot,ptop) ; kg/(m*s) 
VQ_vint = 1./g * 100 * vibeta (interp_levels,VQ_PL(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,ps,pbot,ptop)
copy_VarCoords(ps,UQ_vint)
copy_VarCoords(ps,VQ_vint)
UQ_vint@units = "kg m-1 s-1"
VQ_vint@units = "kg m-1 s-1"

; inter from xy to ll
; New lat lon 
lat = fspan(25,41,33);
lon = fspan(74,105,63);

Q_PL_ll  		= rcm2rgrid(lat2d,lon2d,Q_PL,lat,lon,0)
W_PL_ll  		= rcm2rgrid(lat2d,lon2d,W_PL,lat,lon,0)
Umet_PL_ll  	= rcm2rgrid(lat2d,lon2d,Umet_PL,lat,lon,0)
Vmet_PL_ll  	= rcm2rgrid(lat2d,lon2d,Vmet_PL,lat,lon,0)
height_PL_ll  	= rcm2rgrid(lat2d,lon2d,height_PL,lat,lon,0)
UQ_PL_ll  		= rcm2rgrid(lat2d,lon2d,UQ_PL,lat,lon,0)
VQ_PL_ll  		= rcm2rgrid(lat2d,lon2d,VQ_PL,lat,lon,0)
UQ_vint_ll  	= rcm2rgrid(lat2d,lon2d,UQ_vint,lat,lon,0)
VQ_vint_ll  	= rcm2rgrid(lat2d,lon2d,VQ_vint,lat,lon,0)
pw_ll  			= rcm2rgrid(lat2d,lon2d,pw,lat,lon,0)
ps_ll  			= rcm2rgrid(lat2d,lon2d,ps,lat,lon,0)

copy_VarMeta(Q_PL, Q_PL_ll)
copy_VarMeta(W_PL, W_PL_ll)
copy_VarMeta(Umet_PL, Umet_PL_ll)
copy_VarMeta(Vmet_PL, Vmet_PL_ll)
copy_VarMeta(height_PL, height_PL_ll)
copy_VarMeta(UQ_PL, UQ_PL_ll)
copy_VarMeta(VQ_PL, VQ_PL_ll)
copy_VarMeta(UQ_vint, UQ_vint_ll)
copy_VarMeta(VQ_vint, VQ_vint_ll)
copy_VarMeta(pw, pw_ll)
copy_VarMeta(ps, ps_ll)

; Cal divergenve by UQ VQ
Qdiv_ll = uv2dv_cfd(UQ_PL_ll, VQ_PL_ll, lat, lon,3) 
Qdiv_ll@units = "kg/(kg-s)"
copy_VarCoords(UQ_PL_ll, Qdiv_ll)

; out time name
t_name_out = t_name_in

; output ll
system("rm -rf "+outpath + "wrfpress_daily_ll_d02_" + t_name_out)
outfile = addfile(outpath + "wrfpress_daily_ll_d02_" + t_name_out , "c")

outfile->Times=Times
outfile->lat=lat
outfile->lon=lon
outfile->interp_levels=interp_levels
outfile->Q_PL_ll=Q_PL_ll
outfile->W_PL_ll=W_PL_ll
outfile->Umet_PL_ll=Umet_PL_ll
outfile->Vmet_PL_ll=Vmet_PL_ll
outfile->height_PL_ll=height_PL_ll
outfile->UQ_PL_ll=UQ_PL_ll
outfile->VQ_PL_ll=VQ_PL_ll
outfile->UQ_vint_ll=UQ_vint_ll
outfile->VQ_vint_ll=VQ_vint_ll
outfile->Qdiv_ll=Qdiv_ll
outfile->pw_ll=pw_ll
outfile->ps_ll=ps_ll

; output xy
system("rm -rf "+outpath + "wrfpress_daily_xy_d02_" + t_name_out)
outfile = addfile(outpath + "wrfpress_daily_xy_d02_" + t_name_out , "c")
outfile->Umet_PL=Umet_PL
outfile->Vmet_PL=Vmet_PL
outfile->height_PL=height_PL
outfile->UQ_vint=UQ_vint
outfile->VQ_vint=VQ_vint
outfile->pw=pw


end do







